Sequence = input('Gimme a sample alignment:','s');
fid = fopen(Sequence);
PreAlignment = textscan(fid, '%s');
fclose(fid);
AlignmentACGT = char(PreAlignment{1}); %formatting kludge.
[Maxx, Width] = size(AlignmentACGT);
%  Maxx is the # of sites and Width is the site size.
%  Read in the alignment and convert from acgt-space to XYZ-space.
%  Also convert the acgt-sequence info into its corresponding "binary 
%  matrix" format.
A = [1 -1 -1]; C = [-1 1 -1]; G = [-1 -1 1]; T = [1 1 1]; %see Gube below!
Abine = [1 0 0 0]; Cbine = [0 1 0 0]; Gbine = [0 0 1 0]; Tbine = [0 0 0 1];
AlignmentXYZ = []; BineAlignACGT = [];
for i = (1:Maxx) %down the alignment.
   Tempo1 = []; Tempo2 = [];
   for j = (1:Width) %across each site.
      x = AlignmentACGT(i,j);
      if strcmpi('a', x)
         Tempo1 = [Tempo1 A];
         Tempo2 = [Tempo2 Abine];
      elseif strcmpi('c', x)
         Tempo1 = [Tempo1 C];
         Tempo2 = [Tempo2 Cbine];
      elseif strcmpi('g', x)
         Tempo1 = [Tempo1 G];
         Tempo2 = [Tempo2 Gbine];
      elseif strcmpi('t', x)
         Tempo1 = [Tempo1 T];
         Tempo2 = [Tempo2 Tbine];
      else error('unknown base type in the input file!')
      end
   end
   AlignmentXYZ = [AlignmentXYZ; Tempo1];
   BineAlignACGT = [BineAlignACGT; Tempo2];
end
H = eye(Width*3);
AL = -AlignmentXYZ;
b = -ones(Maxx,1);
options = optimset('LargeScale','off','MaxIter',1e15);
PossumXYZ = quadprog(H,[],AL,b,[],[],[],[],[],options); %in XYZ-space!
Gube = [A; C; G; T]; %this matrix defines Gary's XYZ-space.
BlkGube = []; %transform from XYZ- to acgt-space.
for j = (1:Width)
   BlkGube = blkdiag(BlkGube,Gube);
end
PossumACGT = reshape(BlkGube*PossumXYZ,4,Width); %THE ANSWER!
ScoresACGTpossum = BineAlignACGT*PossumACGT(:);
%mean(PossumACGT) should be 0-vector.
%min(ScoresACGTpossum) should be 1.0.
%
%Print to screen PossumACGT with column and row headers.
JustForDisplay = num2str(vertcat(1:Width,PossumACGT),3);
JustForDisplay = horzcat(['  ';'A ';'C ';'G ';'T '],JustForDisplay);
disp('This is the binding matrix PossumACGT:');
disp(JustForDisplay);
